library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Warning: package 'dplyr' was built under R version 3.4.2
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag(): dplyr, stats
library(ggplot2)
library(LearnBayes)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
Una compañía farmacéutica afirma que su nueva medicina incrementa la probabilidad de concebir un niño (sexo masculino), pero aún no publican estudios. Supón que conduces un experimento en el cual 50 parejas se seleccionan de manera aleatoria de la población, toman la medicina y conciben un bebé, nacen 30 niños y 20 niñas.
a=30
b=20
theta<- rbeta(1000,30,20)
prior<-dbeta(theta,30,20)
base <- ggplot(data_frame(x = c(0, 1)), aes(x))
base +
stat_function(fun = dbeta, args = list(shape1 = 30, shape2 = 20),
aes(colour = "Distribucion inicial"))
N <- 50 # casos
z <- 30 # exitos
final<- dbeta(theta, z+a,N-z+b)/dbeta(theta,a,b)
# Verosimilitud
Like <- theta ^ z * (1 - theta) ^ (N - z)
product <- Like * prior
post <- product / sum(product)
p1 <- base +
stat_function(fun = dbeta, args = list(shape1 = a, shape2 = b),
aes(colour = "inicial"), show.legend = TRUE) +
stat_function(fun = dbeta, args = list(shape1 = z + 1, shape2 = N - z + 1),
aes(colour = "verosimilitud"), show.legend = TRUE) +
stat_function(fun = dbeta, args = list(shape1 = a + z, shape2 = N - z + b),
aes(colour = "posterior"), show.legend = TRUE) +
labs(y = "", colour = "", x = expression(theta))
p1
prior_1<- dbeta(theta,60,40)
prior_2<-dbeta(theta,50,50)
product_1 <- Like * prior_1
product_2 <- Like * prior_2
post_1 <- product_1 / sum(product_1)
post_2 <- product_2 / sum(product_2)
#plot(theta,post_1)
#plot(theta,post_2)
a=60
b=40
p1 <- base +
stat_function(fun = dbeta, args = list(shape1 = a, shape2 = b),
aes(colour = "inicial"), show.legend = FALSE) +
stat_function(fun = dbeta, args = list(shape1 = z + 1, shape2 = N - z + 1),
aes(colour = "verosimilitud"), show.legend = FALSE) +
stat_function(fun = dbeta, args = list(shape1 = a + z, shape2 = N - z + b),
aes(colour = "posterior"), show.legend = FALSE) +
labs(y = "", colour = "", x = expression(theta)) + ggtitle("Beta inicial (60,40)")
a = 50; b = 50
p2 <- base +
stat_function(fun = dbeta, args = list(shape1 = a, shape2 = b),
aes(colour = "inicial")) +
stat_function(fun = dbeta, args = list(shape1 = z + 1, shape2 = N - z + 1),
aes(colour = "verosimilitud")) +
stat_function(fun = dbeta, args = list(shape1 = a + z, shape2 = N - z + b),
aes(colour = "posterior")) +
labs(y = "", colour = "", x = expression(theta)) + ggtitle("Beta inicial (50,50)")
grid.arrange(p1, p2, nrow = 1, widths = c(0.38, 0.62))
Calcula la distribución posterior p(θ|x)∝p(x|θ)p(θ), usando la inicial y verosimilitud que definimos arriba. Una vez que realices la multiplicación debes identificar el núcleo de una distribución Normal, ¿cuáles son sus parámetros (media y varianza)?
2_1
2_2
En el ejercicio anterior hiciste cálculos para el caso de una sola observación. En este ejercicio consideramos el caso en que observamos una muestra x={x1,…,xN}.Crea una función prior que reciba los parámetros μ y τ que definen tus creencias del parámetro desconocido θ y devuelva p(θ), donde p(θ) tiene distribución N(μ,σ2)
prior <- function(mu, tau){
function(theta){
p_0<- (1/sqrt(2*pi*(tau^2)))*exp(-((theta-mu)^2)/(2*(tau^2)))
return(p_0)
}
}
mi_prior<-prior(150,15)
likeNorm <- function(sigma, S, S2, N){
function(theta){
like<- (1/(2*pi*(sigma^(N/2))))*(-(1/(2*sigma^2))*(S2-2*theta*S+N*theta^2))
return(like)
}
}
mi_like<-likeNorm(20,13000,1700000,100)
postRelProb <- function(theta){
mi_like(theta) * mi_prior(theta)
}
# Datos observados
# para cada paso decidimos el movimiento de acuerdo a la siguiente función
caminaAleat <- function(theta){ # theta: valor actual
salto_prop <- rnorm(1, 0, 5) # salto propuesto
theta_prop <- theta + salto_prop # theta propuesta
u <- runif(1)
p_move <- min(postRelProb(theta_prop) / postRelProb(theta), 1) # prob mover
if(p_move>u){
return(theta_prop) # aceptar valor propuesto
}
else{
return(theta) # rechazar
}
}
set.seed(47405)
pasos <- 6000
camino <- numeric(pasos) # vector que guardará las simulaciones
camino[1] <- 15 # valor inicial
# Generamos la caminata aleatoria
for (j in 2:pasos){
camino[j] <- caminaAleat(camino[j - 1])
}
caminata <- data.frame(pasos = 1:pasos, theta = camino)
ggplot(caminata[1:3000, ], aes(x = pasos, y = theta)) +
geom_point(size = 0.8) +
geom_path(alpha = 0.5) +
scale_y_continuous(expression(theta)) +
scale_x_continuous("Tiempo")
camino <- camino[-(1:67)]
post <- lapply(camino,postRelProb) %>% flatten_dbl()
hist(post,breaks=100)